#set working directory
#setwd("~/file_path")

#TERGM -------------------------------------------------------------------

library(missForest)
library(statnet)
library(geosphere)
library(Hmisc)

load("amen_brother/events.RData")
ab_events <- events
load("amen_brother/artist_info.RData")
ab_artist_info <- artist_info
load("amen_brother/spotify_data.RData")
ab_spotify_data <- spotify_data
load("amen_brother/geo_thesaurus.RData")
ab_geo_thesaurus <- geo_thesaurus
load("amen_brother/discogs_ids.RData")
ab_discogs_ids <- discogs_ids
load("amen_brother/years.RData")
ab_years <- years
rm(list = c("events", "artist_info", "spotify_data", "geo_thesaurus", "discogs_ids", "years"))

load("think_about_it/events.RData")
tai_events <- events
load("think_about_it/artist_info.RData")
tai_artist_info <- artist_info
load("think_about_it/spotify_data.RData")
tai_spotify_data <- spotify_data
load("think_about_it/geo_thesaurus.RData")
tai_geo_thesaurus <- geo_thesaurus
load("think_about_it/discogs_ids.RData")
tai_discogs_ids <- discogs_ids
load("think_about_it/years.RData")
tai_years <- years
rm(list = c("events", "artist_info", "spotify_data", "geo_thesaurus", "discogs_ids", "years"))

load("funky_drummer/events.RData")
fd_events <- events
load("funky_drummer/artist_info.RData")
fd_artist_info <- artist_info
load("funky_drummer/spotify_data.RData")
fd_spotify_data <- spotify_data
load("funky_drummer/geo_thesaurus.RData")
fd_geo_thesaurus <- geo_thesaurus
load("funky_drummer/discogs_ids.RData")
fd_discogs_ids <- discogs_ids
load("funky_drummer/years.RData")
fd_years <- years
rm(list = c("events", "artist_info", "spotify_data", "geo_thesaurus", "discogs_ids", "years"))

events <- rbind(ab_events, tai_events, fd_events)
events <- events[order(events$years), ]

artist_info <- rbind(ab_artist_info, tai_artist_info, fd_artist_info)
spotify_data <- rbind(ab_spotify_data, tai_spotify_data, fd_spotify_data)
geo_thesaurus <- rbind(ab_geo_thesaurus, tai_geo_thesaurus, fd_geo_thesaurus)
discogs_ids <- c(ab_discogs_ids, tai_discogs_ids, fd_discogs_ids)
years <- c(ab_years, tai_years, fd_years)

#remove everything besides needed objects
rm(list = ls()[-which(ls() %in% c("events", "artist_info", "spotify_data", "geo_thesaurus", "discogs_ids", "years"))])

verts <- sort(as.numeric(unique(c(as.character(events[which(events$years %in% 1984:2017), 1]), as.character(events[which(events$years %in% 1984:2017), 2])))))

#generate matrices for each year
i <- 1984
while(i <= 2017){
  #generate empty labelled matrix
  temp_matrix <- matrix(0, nrow = length(verts), ncol = length(verts))
  rownames(temp_matrix) <- verts
  colnames(temp_matrix) <- verts
  
  #store events from target year
  temp_events <- events[which(events$years == i), 1:2]
  
  #add events to matrix
  j <- 1
  while(j <= nrow(temp_events)){
    temp <- as.numeric(as.character(unlist(temp_events[j, ])))
    temp_matrix[which(rownames(temp_matrix) == temp[1]), which(rownames(temp_matrix) == temp[2])] <- temp_matrix[which(rownames(temp_matrix) == temp[1]), which(rownames(temp_matrix) == temp[2])] + 1
    temp_matrix[which(rownames(temp_matrix) == temp[2]), which(rownames(temp_matrix) == temp[1])] <- temp_matrix[which(rownames(temp_matrix) == temp[2]), which(rownames(temp_matrix) == temp[1])] + 1
    j <- j + 1
  }
  
  temp_network <- as.network(temp_matrix, matrix.type = "adjacency", directed = FALSE, loops = FALSE)
  temp_network %e% "weight" <- temp_matrix
  
  #store and advance i
  assign(paste("net", i, sep = ""), temp_network)
  i <- i + 1
}

nets <- list(net1984, net1985, net1986, net1987, net1988, 
             net1989, net1990, net1991, net1992, net1993, 
             net1994, net1995, net1996, net1997, net1998,
             net1999, net2000, net2001, net2002, net2003,
             net2004, net2005, net2006, net2007, net2008,
             net2009, net2010, net2011, net2012, net2013,
             net2014, net2015, net2016, net2017)

gender <- c()
popularity <- c()
followers <- c()
i <- 1
while(i <= length(verts)){
  gender <- c(gender, artist_info$gender[which(artist_info$discogs == verts[i])][1])
  popularity <- c(popularity, spotify_data$popularity[which(artist_info$discogs == verts[i])][1])
  followers <- c(followers, spotify_data$followers[which(artist_info$discogs == verts[i])][1])
  
  i <- i + 1
}

#female = -1; male = 1; 0 = NA or other
gender[which(gender == 1)] <- -1
gender[which(gender == 2)] <- 1
gender[which(gender == 3)] <- 0
gender[which(gender == 4)] <- 0
gender[is.na(gender)] <- 0

#retrieve lon and lat for each artist
lon <- c()
lat <- c()
i <- 1
while(i <= length(verts)){
  lon <- c(lon, geo_thesaurus$lon[which(geo_thesaurus$id == verts[i])][1])
  lat <- c(lat, geo_thesaurus$lat[which(geo_thesaurus$id == verts[i])][1])
  i <- i + 1
}

coords <- matrix(c(lon, lat), nrow = length(lon), ncol = 2)

#calculate mean distance of each point from all other points (in meters) on an ellipsoid, excluding individuals without geographic info...
i <- 1
meandist <- c()
while(i <= length(lon)){
  if(!is.na(lon[i])){
    temp <- coords[which(!is.na(lon[-1])) + 1, ]
    meandist <- c(meandist, mean(distGeo(coords[i,], temp)))
  }
  else{
    meandist <- c(meandist, NA)
  }
  i <- i + 1
}

vars <- data.frame(popularity, followers, meandist)
vars_imputed <- missForest(vars)

popularity <- vars_imputed$ximp$popularity
followers <- vars_imputed$ximp$followers
meandist <- vars_imputed$ximp$meandist

#center continuous variables around 0
popularity <- c(scale(popularity, center = TRUE, scale = FALSE))
followers <- c(scale(followers, center = TRUE, scale = FALSE))
meandist <- c(scale(meandist, center = TRUE, scale = FALSE))

i <- 1
while(i <= length(nets)){
  set.vertex.attribute(nets[[i]], "meandist", meandist)
  set.vertex.attribute(nets[[i]], "popularity", popularity)
  set.vertex.attribute(nets[[i]], "followers", followers)
  set.vertex.attribute(nets[[i]], "gender", gender)
  
  i <- i + 1
}

#compare all possible models for pre-2000 and post-2000
stergm_1_pre <- stergm(nets, formation = ~edges, dissolution = ~edges, estimate = "CMLE", times = c(1:16), verbose = FALSE)
stergm_2_pre <- stergm(nets, formation = ~edges+nodecov("meandist"), dissolution = ~edges, estimate = "CMLE", times = c(1:16), verbose = FALSE)
stergm_3_pre <- stergm(nets, formation = ~edges+absdiff("popularity"), dissolution = ~edges, estimate = "CMLE", times = c(1:16), verbose = FALSE)
stergm_4_pre <- stergm(nets, formation = ~edges+absdiff("followers"), dissolution = ~edges, estimate = "CMLE", times = c(1:16), verbose = FALSE)
stergm_5_pre <- stergm(nets, formation = ~edges+nodematch("gender", diff = TRUE), dissolution = ~edges, estimate = "CMLE", times = c(1:16), verbose = FALSE)
stergm_6_pre <- stergm(nets, formation = ~edges+nodecov("meandist")+absdiff("popularity"), dissolution = ~edges, estimate = "CMLE", times = c(1:16), verbose = FALSE)
stergm_7_pre <- stergm(nets, formation = ~edges+nodecov("meandist")+absdiff("followers"), dissolution = ~edges, estimate = "CMLE", times = c(1:16), verbose = FALSE)
stergm_8_pre <- stergm(nets, formation = ~edges+nodecov("meandist")+nodematch("gender", diff = TRUE), dissolution = ~edges, estimate = "CMLE", times = c(1:16), verbose = FALSE)
stergm_9_pre <- stergm(nets, formation = ~edges+absdiff("popularity")+absdiff("followers"), dissolution = ~edges, estimate = "CMLE", times = c(1:16), verbose = FALSE)
stergm_10_pre <- stergm(nets, formation = ~edges+absdiff("popularity")+nodematch("gender", diff = TRUE), dissolution = ~edges, estimate = "CMLE", times = c(1:16), verbose = FALSE)
stergm_11_pre <- stergm(nets, formation = ~edges+absdiff("followers")+nodematch("gender", diff = TRUE), dissolution = ~edges, estimate = "CMLE", times = c(1:16), verbose = FALSE)
stergm_12_pre <- stergm(nets, formation = ~edges+nodecov("meandist")+absdiff("popularity")+absdiff("followers"), dissolution = ~edges, estimate = "CMLE", times = c(1:16), verbose = FALSE)
stergm_13_pre <- stergm(nets, formation = ~edges+nodecov("meandist")+absdiff("popularity")+nodematch("gender", diff = TRUE), dissolution = ~edges, estimate = "CMLE", times = c(1:16), verbose = FALSE)
stergm_14_pre <- stergm(nets, formation = ~edges+nodecov("meandist")+absdiff("followers")+nodematch("gender", diff = TRUE), dissolution = ~edges, estimate = "CMLE", times = c(1:16), verbose = FALSE)
stergm_15_pre <- stergm(nets, formation = ~edges+absdiff("popularity")+absdiff("followers")+nodematch("gender", diff = TRUE), dissolution = ~edges, estimate = "CMLE", times = c(1:16), verbose = FALSE)
stergm_16_pre <- stergm(nets, formation = ~edges+nodecov("meandist")+absdiff("popularity")+absdiff("followers")+nodematch("gender", diff = TRUE), dissolution = ~edges, estimate = "CMLE", times = c(1:16), verbose = FALSE)

#get AIC values...
i <- 1
aic_pre <- c()
while(i <= 16){
  aic_pre <- c(aic_pre, min(summary(get(print(paste("stergm_", i, "_pre", sep = ""))))$formation$aic))
  i <- i + 1
}

#find best models (deltaAIC < 2)...
delta_aic_pre <- aic_pre-min(aic_pre)
best_pre <- which(delta_aic_pre < 2)

#output AIC table
var_ref <- c("nodematch.gender.0", "absdiff.popularity", "absdiff.followers", "nodecov.meandist")
ilv <- list()
i <- 1
while(i <= 16){
  temp <- row.names(summary(get(paste0("stergm_", i, "_pre")))$formation$asycov)
  temp <- paste(as.character(which(var_ref %in% temp)), collapse = " ")
  if(nchar(temp) == 0){
    temp <- c("0")
  }
  ilv[i] <- temp
  i <- i + 1
}

stergm_table_pre <- data.frame(unlist(ilv), aic_pre, delta_aic_pre)[order(delta_aic_pre),]
colnames(stergm_table_pre) <- c("ILVs", "AIC", "deltaAIC")

#goodness-of-fit analysis
gof_pre <- gof(stergm_16_pre)

stergm_1_post <- stergm(nets, formation = ~edges, dissolution = ~edges, estimate = "CMLE", times = c(17:34), verbose = FALSE)
stergm_2_post <- stergm(nets, formation = ~edges+nodecov("meandist"), dissolution = ~edges, estimate = "CMLE", times = c(17:34), verbose = FALSE)
stergm_3_post <- stergm(nets, formation = ~edges+absdiff("popularity"), dissolution = ~edges, estimate = "CMLE", times = c(17:34), verbose = FALSE)
stergm_4_post <- stergm(nets, formation = ~edges+absdiff("followers"), dissolution = ~edges, estimate = "CMLE", times = c(17:34), verbose = FALSE)
stergm_5_post <- stergm(nets, formation = ~edges+nodematch("gender", diff = TRUE), dissolution = ~edges, estimate = "CMLE", times = c(17:34), verbose = FALSE)
stergm_6_post <- stergm(nets, formation = ~edges+nodecov("meandist")+absdiff("popularity"), dissolution = ~edges, estimate = "CMLE", times = c(17:34), verbose = FALSE)
stergm_7_post <- stergm(nets, formation = ~edges+nodecov("meandist")+absdiff("followers"), dissolution = ~edges, estimate = "CMLE", times = c(17:34), verbose = FALSE)
stergm_8_post <- stergm(nets, formation = ~edges+nodecov("meandist")+nodematch("gender", diff = TRUE), dissolution = ~edges, estimate = "CMLE", times = c(17:34), verbose = FALSE)
stergm_9_post <- stergm(nets, formation = ~edges+absdiff("popularity")+absdiff("followers"), dissolution = ~edges, estimate = "CMLE", times = c(17:34), verbose = FALSE)
stergm_10_post <- stergm(nets, formation = ~edges+absdiff("popularity")+nodematch("gender", diff = TRUE), dissolution = ~edges, estimate = "CMLE", times = c(17:34), verbose = FALSE)
stergm_11_post <- stergm(nets, formation = ~edges+absdiff("followers")+nodematch("gender", diff = TRUE), dissolution = ~edges, estimate = "CMLE", times = c(17:34), verbose = FALSE)
stergm_12_post <- stergm(nets, formation = ~edges+nodecov("meandist")+absdiff("popularity")+absdiff("followers"), dissolution = ~edges, estimate = "CMLE", times = c(17:34), verbose = FALSE)
stergm_13_post <- stergm(nets, formation = ~edges+nodecov("meandist")+absdiff("popularity")+nodematch("gender", diff = TRUE), dissolution = ~edges, estimate = "CMLE", times = c(17:34), verbose = FALSE)
stergm_14_post <- stergm(nets, formation = ~edges+nodecov("meandist")+absdiff("followers")+nodematch("gender", diff = TRUE), dissolution = ~edges, estimate = "CMLE", times = c(17:34), verbose = FALSE)
stergm_15_post <- stergm(nets, formation = ~edges+absdiff("popularity")+absdiff("followers")+nodematch("gender", diff = TRUE), dissolution = ~edges, estimate = "CMLE", times = c(17:34), verbose = FALSE)
stergm_16_post <- stergm(nets, formation = ~edges+nodecov("meandist")+absdiff("popularity")+absdiff("followers")+nodematch("gender", diff = TRUE), dissolution = ~edges, estimate = "CMLE", times = c(17:34), verbose = FALSE)

#get AIC values...
i <- 1
aic_post <- c()
while(i <= 16){
  aic_post <- c(aic_post, min(summary(get(print(paste("stergm_", i, "_post", sep = ""))))$formation$aic))
  i <- i + 1
}

#find best models (deltaAIC < 2)...
delta_aic_post <- aic_post-min(aic_post)
best_post <- which(delta_aic_post < 2)

#output AIC table
var_ref <- c("nodematch.gender.0", "absdiff.popularity", "absdiff.followers", "nodecov.meandist")
ilv <- list()
i <- 1
while(i <= 16){
  temp <- row.names(summary(get(paste0("stergm_", i, "_post")))$formation$asycov)
  temp <- paste(as.character(which(var_ref %in% temp)), collapse = " ")
  if(nchar(temp) == 0){
    temp <- c("0")
  }
  ilv[i] <- temp
  i <- i + 1
}

stergm_table_post <- data.frame(unlist(ilv), aic_post, delta_aic_post)[order(delta_aic_post),]
colnames(stergm_table_post) <- c("ILVs", "AIC", "deltaAIC")

#goodness-of-fit analysis
gof_post <- gof(stergm_16_post)

stergm_pre_1994 <- stergm(nets, formation = ~edges+nodecov("meandist")+absdiff("popularity")+absdiff("followers")+nodematch("gender", diff = TRUE), dissolution = ~edges, estimate = "CMLE", times = c(1:10), verbose = FALSE)
stergm_post_1994 <- stergm(nets, formation = ~edges+nodecov("meandist")+absdiff("popularity")+absdiff("followers")+nodematch("gender", diff = TRUE), dissolution = ~edges, estimate = "CMLE", times = c(11:34), verbose = FALSE)

stergm_pre_1996 <- stergm(nets, formation = ~edges+nodecov("meandist")+absdiff("popularity")+absdiff("followers")+nodematch("gender", diff = TRUE), dissolution = ~edges, estimate = "CMLE", times = c(1:12), verbose = FALSE)
stergm_post_1996 <- stergm(nets, formation = ~edges+nodecov("meandist")+absdiff("popularity")+absdiff("followers")+nodematch("gender", diff = TRUE), dissolution = ~edges, estimate = "CMLE", times = c(13:34), verbose = FALSE)

stergm_pre_1998 <- stergm(nets, formation = ~edges+nodecov("meandist")+absdiff("popularity")+absdiff("followers")+nodematch("gender", diff = TRUE), dissolution = ~edges, estimate = "CMLE", times = c(1:14), verbose = FALSE)
stergm_post_1998 <- stergm(nets, formation = ~edges+nodecov("meandist")+absdiff("popularity")+absdiff("followers")+nodematch("gender", diff = TRUE), dissolution = ~edges, estimate = "CMLE", times = c(15:34), verbose = FALSE)

stergm_pre_2002 <- stergm(nets, formation = ~edges+nodecov("meandist")+absdiff("popularity")+absdiff("followers")+nodematch("gender", diff = TRUE), dissolution = ~edges, estimate = "CMLE", times = c(1:18), verbose = FALSE)
stergm_post_2002 <- stergm(nets, formation = ~edges+nodecov("meandist")+absdiff("popularity")+absdiff("followers")+nodematch("gender", diff = TRUE), dissolution = ~edges, estimate = "CMLE", times = c(19:34), verbose = FALSE)

stergm_pre_2004 <- stergm(nets, formation = ~edges+nodecov("meandist")+absdiff("popularity")+absdiff("followers")+nodematch("gender", diff = TRUE), dissolution = ~edges, estimate = "CMLE", times = c(1:20), verbose = FALSE)
stergm_post_2004 <- stergm(nets, formation = ~edges+nodecov("meandist")+absdiff("popularity")+absdiff("followers")+nodematch("gender", diff = TRUE), dissolution = ~edges, estimate = "CMLE", times = c(21:34), verbose = FALSE)

stergm_pre_2006 <- stergm(nets, formation = ~edges+nodecov("meandist")+absdiff("popularity")+absdiff("followers")+nodematch("gender", diff = TRUE), dissolution = ~edges, estimate = "CMLE", times = c(1:22), verbose = FALSE)
stergm_post_2006 <- stergm(nets, formation = ~edges+nodecov("meandist")+absdiff("popularity")+absdiff("followers")+nodematch("gender", diff = TRUE), dissolution = ~edges, estimate = "CMLE", times = c(23:34), verbose = FALSE)

#number of unique individuals involved in collaborations, not the number of collaborations...
pre_sample_sizes <- c(length(unique(c(events$ind1[which(as.character(events$years) < 1994)], events$ind2[which(as.character(events$years) < 1994)]))),
                      length(unique(c(events$ind1[which(as.character(events$years) < 1996)], events$ind2[which(as.character(events$years) < 1996)]))),
                      length(unique(c(events$ind1[which(as.character(events$years) < 1998)], events$ind2[which(as.character(events$years) < 1998)]))),
                      length(unique(c(events$ind1[which(as.character(events$years) < 2000)], events$ind2[which(as.character(events$years) < 2000)]))),
                      length(unique(c(events$ind1[which(as.character(events$years) < 2002)], events$ind2[which(as.character(events$years) < 2002)]))),
                      length(unique(c(events$ind1[which(as.character(events$years) < 2004)], events$ind2[which(as.character(events$years) < 2004)]))),
                      length(unique(c(events$ind1[which(as.character(events$years) < 2006)], events$ind2[which(as.character(events$years) < 2006)]))))

post_sample_sizes <- c(length(unique(c(events$ind1[which(as.character(events$years) >= 1994)], events$ind2[which(as.character(events$years) >= 1994)]))),
                       length(unique(c(events$ind1[which(as.character(events$years) >= 1996)], events$ind2[which(as.character(events$years) >= 1996)]))),
                       length(unique(c(events$ind1[which(as.character(events$years) >= 1998)], events$ind2[which(as.character(events$years) >= 1998)]))),
                       length(unique(c(events$ind1[which(as.character(events$years) >= 2000)], events$ind2[which(as.character(events$years) >= 2000)]))),
                       length(unique(c(events$ind1[which(as.character(events$years) >= 2002)], events$ind2[which(as.character(events$years) >= 2002)]))),
                       length(unique(c(events$ind1[which(as.character(events$years) >= 2004)], events$ind2[which(as.character(events$years) >= 2004)]))),
                       length(unique(c(events$ind1[which(as.character(events$years) >= 2006)], events$ind2[which(as.character(events$years) >= 2006)]))))

pre_table <- data.frame(summary(stergm_pre_1994)$formation$coefs[,1], summary(stergm_pre_1994)$formation$coefs[,4],
                        summary(stergm_pre_1996)$formation$coefs[,1], summary(stergm_pre_1996)$formation$coefs[,4],
                        summary(stergm_pre_1998)$formation$coefs[,1], summary(stergm_pre_1998)$formation$coefs[,4],
                        summary(stergm_16_pre)$formation$coefs[,1], summary(stergm_16_pre)$formation$coefs[,4],
                        summary(stergm_pre_2002)$formation$coefs[,1], summary(stergm_pre_2002)$formation$coefs[,4],
                        summary(stergm_pre_2004)$formation$coefs[,1], summary(stergm_pre_2004)$formation$coefs[,4],
                        summary(stergm_pre_2006)$formation$coefs[,1], summary(stergm_pre_2006)$formation$coefs[,4])
                        
pre_table <- pre_table[c(-1, -6),]
pre_table <- format(pre_table, digits = 2, scientific = TRUE)
pre_table <- rbind(c("Estimate", "p-value", "Estimate", "p-value", "Estimate", "p-value", "Estimate", "p-value", "Estimate", "p-value", "Estimate", "p-value", "Estimate", "p-value"), pre_table)
pre_table <- rbind(c(paste0("n = ", pre_sample_sizes[1]), "", paste0("n = ", pre_sample_sizes[2]), "", paste0("n = ", pre_sample_sizes[3]), "", paste0("n = ", pre_sample_sizes[4]), "", paste0("n = ", pre_sample_sizes[5]), "", paste0("n = ", pre_sample_sizes[6]), "", paste0("n = ", pre_sample_sizes[7]), ""), pre_table)
colnames(pre_table) <- c("Pre-1994", "", "Pre-1996", "", "Pre-1998", "", "Pre-2000", "", "Pre-2002", "", "Pre-2004", "", "Pre-2006", "")
rownames(pre_table) <- c("", "_", "Mean Distance", "Popularity", "Followers", "Gender (F)", "Gender (M)")

post_table <- data.frame(summary(stergm_post_1994)$formation$coefs[,1], summary(stergm_post_1994)$formation$coefs[,4],
                        summary(stergm_post_1996)$formation$coefs[,1], summary(stergm_post_1996)$formation$coefs[,4],
                        summary(stergm_post_1998)$formation$coefs[,1], summary(stergm_post_1998)$formation$coefs[,4],
                        summary(stergm_16_post)$formation$coefs[,1], summary(stergm_16_post)$formation$coefs[,4],
                        summary(stergm_post_2002)$formation$coefs[,1], summary(stergm_post_2002)$formation$coefs[,4],
                        summary(stergm_post_2004)$formation$coefs[,1], summary(stergm_post_2004)$formation$coefs[,4],
                        summary(stergm_post_2006)$formation$coefs[,1], summary(stergm_post_2006)$formation$coefs[,4])

post_table <- post_table[c(-1, -6),]
post_table <- format(post_table, digits = 2, scientific = TRUE)
post_table <- rbind(c("Estimate", "p-value", "Estimate", "p-value", "Estimate", "p-value", "Estimate", "p-value", "Estimate", "p-value", "Estimate", "p-value", "Estimate", "p-value"), post_table)
post_table <- rbind(c(paste0("n = ", post_sample_sizes[1]), "", paste0("n = ", post_sample_sizes[2]), "", paste0("n = ", post_sample_sizes[3]), "", paste0("n = ", post_sample_sizes[4]), "", paste0("n = ", post_sample_sizes[5]), "", paste0("n = ", post_sample_sizes[6]), "", paste0("n = ", post_sample_sizes[7]), ""), post_table)
colnames(post_table) <- c("Post-1994", "", "Post-1996", "", "Post-1998", "", "Post-2000", "", "Post-2002", "", "Post-2004", "", "Post-2006", "")
rownames(post_table) <- c("", "_", "Mean Distance", "Popularity", "Followers", "Gender (F)", "Gender (M)")
